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Abstract 

Numerical  simulations  using  the  lattice  Boltzmann  method  are  presented 
for  the  two-  and  three-dimensional  decaying  homogeneous  isotropic  turbu- 
lence for  low,  medium  and  high  Reynolds  numbers.  Time  history  of  global 
statistical  quantities,  wave  number  spectra,  and  vorticity  contour  plots  are 
compared  with  those  of  the  higher-order  method  of  lines.  Comparisons  be- 
tween the  square  lattice  and  the  triangular  (FHP)  lattice  models  are  also 
performed.  It  is  found  that  the  lattice  Boltzmann  method  is  able  to  re- 
produce the  dynamics  of  decaying  turbulence  and  could  be  an  alternative 
for  solving  the  Navier-Stokes  equations.  Computational  costs  of  the  lattice 
Boltzmann  method  is  less  than  half  of  that  of  the  lOth-order  method  of 
lines. 

1.  Introduction 

The  rapid  development  and  introduction  of  new  supercomputer  systems 
over  the  last  decade  has  opened  new  opportunities  for  numerical  studies  of 
incompressible  fluid  flows.  The  direct  numerical  simulations  of  turbulence  is 
one  of  such  problems.  So  far  almost  all  direct  simulations  of  turbulence  has 
been  carried  out  by  either  spectral[l]  or  pseudo  spectral[2,3]  approximation 
to  spatial  derivatives.  However,  these  methods  which  require  the  use  of 
series  are  global  in  character  so  that  they  are  quite  unsuitable  for  complex 
geometry  problems  and  for  parallel  computing.  Therefore,  the  development 
of  more  flexible  and  efficient  methods  is  hoped  for  in  the  simulations  of 
turbulence. 
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We  have  proposed  a new  higher-order  method  based  on  a method  of 
lines  (MOL)  approach[4-7]  and  demonstrated  that  results  obtained  by  the 
method  were  comparable  to  those  using  the  pseudospectral  method  with 
less  than  one  sixth  of  the  computational  time  in  direct  simulation  of  two- 
dimensional  homogeneous  isotropic  turbulence  on  513  x 513  grid  points. 

In  the  later  half  of  the  80’s,  a novel  technique  called  Lattice  Gas  Au- 
tomata (LGA)  for  solving  the  Navier-Stokes  equations  was  developed.  Since 
the  first  two-dimensional  model  representing  incompressible  Navier-Stokes 
equations  was  proposed  by  Frisch,  Hasslacher,  and  Pomeau  (FHP)  in  1986[8], 
LGA  have  attracted  much  attention  as  promising  methods  for  solving  a va- 
riety of  partial  differential  equations  and  modeling  physical  phenomena. 
The  basic  idea  of  LGA  methods  is  to  represent  the  fluid  as  an  ensemble 
of  interacting  low-order  bit-computers  situated  at  regularly  spaced  lattice 
nodes.  In  the  FHP  model,  the  underlying  lattice  is  a close-packed  equi- 
lateral triangular  lattice  with  nodes  at  triangle  vertices,  each  node  has  a 
seven-bit  state  with  the  first  bit  specifying  the  existence  or  not  of  a par- 
ticle at  rest  at  the  lattice  node  and  the  remaining  six  bits  specifying  the 
presence  or  not  of  a particle  traveling  at  an  angle  dj  — (7r/3)y  (0  < j < 6) 
along  the  legs  of  the  triangular  lattice.  Each  particle  (except  a rest  parti- 
cle) moves  one  lattice  distance  in  one  fundamental  time  interval.  After  the 
particles  propagate  they  then  interact  according  to  certain  collision  rules. 
Although  the  LGA  method  has  provided  a fast  and  efficient  way  for  solv- 
ing partial  differential  equations,  there  exist  some  fundamental  problems 
in  this  method  in  simulating  realistic  fluid  flows  obeying  the  Navier-Stokes 
equations.  Besides  its  intrinsic  noisy  character  which  makes  the  computa- 
tional accuracy  difficult  to  achieve,  it  contains  certain  properties  even  in  the 
fluid  limit.  The  lattice  gas  fluid  momentum  equations  cannot  be  reduced 
to  the  Navier-Stokes  equations  because  of  two  fundamental  problems.  The 
first  is  the  non-Galilean  invariance  property  due  to  the  density  dependence 
of  the  convection  coefficient.  This  limits  the  validity  of  the  LGA  method 
only  a strict  incompressible  region.  Second  the  pressure  has  an  explicit  and 
unphysical  velocity  dependence.  To  avoid  some  of  those  inherent  problems, 
several  lattice  Boltzmann  (LB)  models  have  been  proposed[9-15].  The  main 
feature  of  the  LB  method  is  to  replace  the  particle  occupation  variables  rii 
(Boolean  variables)  by  the  single-particle  distribution  functions  (real  vari- 
ables) fi  = (rij),  where  ( ) denotes  a local  ensemble  average,  in  the  evolution 
equation,  i.e.,  the  lattice  Boltzmann  equation.  The  LB  method  as  a numer- 
ical scheme  was  first  proposed  by  McNamara  and  Zaretti[9].  In  their  model, 
the  form  of  collision  operator  is  the  same  as  in  the  LGA,  written  in  terms 
of  distribution  functions  and  completely  neglecting  the  effect  of  correla- 
tions between  particles.  Higuera,  Jimenez,  and  Succi[10,ll]  introduced  a 
linearized  collision  operator  that  is  a matrix  and  has  no  correspondence  to 
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the  detailed  collision  rules.  Statistical  noise  is  completely  eliminated  in  both 
models;  however,  the  other  problems  remain,  since  the  equilibrium  distri- 
bution is  still  Fermi-Diracs.  The  LB  model  proposed  by  Chen  et  al[12,13] 
and  Qian  et  al[14,15]  abandons  Fermi-Dirac  statistic  and  applies  the  single 
relaxation  time  approximation  first  introduced  by  Bhatnager,  Gross,  and 
Krook  in  1954[16],  to  greatly  simplify  the  collision  operator.  This  model  is 
called  the  lattice  BGK  (LBGK)  model. 

This  paper  organized  as  follows.  In  section  2 the  lattice  Boltzmann 
methods  simulating  the  Navier-Stokes  equations  are  discussed.  The  lat- 
tice Boltzmann  simulation  of  two-dimensional  homogeneous  isotropic  tur- 
bulence is  presented  in  section  3.  Three-dimensional  homogeneous  isotropic 
turbulence  is  presented  in  section  4.  Accuracy  and  efficiency  of  the  lattice 
Boltzmann  method  in  comparison  with  the  conventional  higher-order  MOL 
approach  are  also  discussed.  The  final  section  contains  concluding  remarks. 

2.  Lattice  Boltzmann  Method 

2.1.  SQUARE  LATTICE  MODEL 

In  this  section  an  outline  is  given  of  the  LB  methods  with  BGK  model 
for  the  collision  operator.  A square  lattice  with  unit  spacing  is  used  on 
which  each  node  has  eight  nearest  neighbors  connected  by  eight  links  as 
shown  in  Fig.l.  Particles  can  only  reside  on  the  nodes  and  move  to  their 
nearest  neighbors  along  these  links  in  the  unit  time.  Hence,  there  are  two 
types  of  moving  particles.  Particles  of  type  1 move  along  the  axes  with 
speed  |epj|  = 1 and  particle  of  type  2 move  along  the  diagonal  directions 
with  speed  |e2,j|  = V2.  Rest  particles  with  speed  zero  are  also  allowed  at 
each  node.  The  occupation  of  the  three  types  of  particles  is  represented 
by  the  single-particle  distribution  function,  fai(x,t),  where  subscripts  a 
and  i indicate  the  type  of  particle  and  the  velocity  direction,  respectively. 
When  <7  = 0,  there  is  only  /oi-  The  distribution  function,  /aj(x,  t),  is  the 
probability  of  finding  a particle  at  node  x and  time  t with  velocity  effj. 
According  to  Bhatnagar,  Gross,  and  Krook  (BGK),  the  collision  operator 
is  simplified  using  the  single  time  relaxation  approximation.  Hence,  the 
lattice  Boltzmann  BGK  (LBGK)  equation  (in  lattice  unit)  is 

/CTi(x  + eai,t+  1)  - Ui{x,t)  = (x,  t)  - /S}(x,t)]  (1) 

T 

where  /^(x,  t)  is  the  equilibrium  distribution  at  x,f  and  t is  the  single 
relaxation  time  which  controls  the  rate  of  approach  to  equilibrium.  The 
density  per  node,  p , and  the  macroscopic  velocity,  u , are  defined  in  terms 
of  the  particle  distribution  function  by 


(2) 
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A suitable  equilibrium  distribution  can  be  chosen  in  the  following  form  for 
particles  of  each  type 

,(0)  2 2 

/oi  = pot  - gpu 


fii]  = pP  + ^P(eii  • u)  + ip(e u • u)2 


1 2 


AO) 

J2i 


(1-40 -a)  1 , , 1 , , 2 1 2 

P 1 + ^P(e2i  - u)  + &P(e2i  ‘ u)  ~ 24 PU 


The  relaxation  time  is  related  to  the  viscosity  by 

2 r - 1 


v = 


6 


(3) 

(4) 

(5) 

(6) 


where  v is  the  kinematic  viscosity  measured  in  lattice  units.  In  ref.  [17],  Hou 
et  al  used  the  value  of  a = 4/9  and  0 = 1/9. 

The  equilibrium  populations  are  determined  by  assuming  that  they  can 
be  expressed  as  a power  series  in  velocity  and  density  of  the  form: 


AO) 
J oi 


■Affi(p)  + S0-j(p)effj  • U + C<ri(p)(®<n  • u)  + .D(p)<jjU 


(7) 


A Chapman-Enskog  procedure  is  then  applied  to  determine  the  macro- 
scopic behavior  of  this  model.  The  values  of  A„i,  BIJl,  Cai  and  I)nl  are  cho- 
sen so  that  the  macroscopic  behavior  matches  the  Navier-Stokes  equations 
to  as  high  an  order  as  possible.  The  resulting  continuity  and  momentum 
equations  follow. 


dp  dpup 
dt  dxp 


+ 0{e2)  = 0 


(8) 


dua 

p~m 


4 -pup 


dua 

dxp 


dp  d f ( dup  dun  \ \ 

dxa  dxp  \9a:a  dxp  J ) 


+ 0(e2)  + 0(M *)  (9) 


Characteristic  dimensionless  parameters  are  the  Mach  number,  Ma  = 
\[ZU j c where  U is  a characteristic  macroscopic  flow  speed  and  c = 1 in 
lattice  units,  the  Knudsen  number  which  is  proportional  to  e = ct/L  where 
L is  a macroscopic  flow  length,  and  the  Reynolds  number,  Re  = pUL/p. 


2.2.  TRIANGULAR  (FHP)  LATTICE  MODEL 

Another  lattice  model  commonly  used  in  two-dimensional  LB  simulation  is 
a triangular  lattice  (FHP)  model.  There  are  two  types  of  particles  on  each 
node  of  the  FHP  model:  rest  particles  and  moving  particles  with  unit  veloc- 
ity ei  along  six  directions  as  shown  in  Fig.2.  The  equilibrium  distributions 
for  the  FHP  model  are  given  as, 

/o(0)  = do-  pu2  = pa-  pu 2 


(10) 
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fi°]  =d+\p  (ei ' u)  + 2(e* ' u)2  ~ ^u2 
= P~jf?  + ^ P (ei  • u)  + 2(e;  • u)2  - ^u2  (11) 


where  a is  an  adjustable  parameter.  If  the  ratio  of  rest  and  moving  parti- 
cles is  defined  as  A = do/d,  the  pressure  is  determined  by  the  isothermal 
equation  of  state, 


(1  ~a)P 
2 


and  the  speed  of  sound  is 


3 

A + 6 


(13) 


The  viscosity  is  related  to  the  relaxation  time  through  an  equation  of  the 
form 

2 t - 1 


2.3.  CUBIC  LATTICE  MODEL 

A cubic  lattice[17]  with  unit  spacing  is  used  on  which  each  node  has  fourteen 
nearest  neighbors  connected  by  fourteen  links.  Particles  can  only  reside  on 
the  nodes  and  move  to  their  nearest  neighbors  along  these  links  in  the  unit 
time  as  shown  in  Fig. 3.  Hence,  there  are  two  types  of  moving  particles. 
Particles  of  type  1 move  along  the  axes  with  speed  |ei,«|  = 1 and  particle 
of  type  2 move  along  the  links  to  the  corners  with  speed  |e2,»|  = \/3-  Rest 
particles  with  speed  zero  are  also  allowed  at  each  node. 

A suitable  equilibrium  distribution  can  be  chosen  in  the  following  form 
for  particles  of  each  type 

/oi)  =pOi~  ~pn2  (15) 

fif  =PP+  |p(ei<  • u)  + ^p{eu  • u)2  - |rpu2  (16) 

fv  = P~ — — — + Y^p(e2z  • u)  + ^p{e2i  ■ u)2  - ^pu2  (17) 

Values  of  a = 2/9  and  (3  — 1/9  are  used.  The  relaxation  time  is  related  to 
the  viscosity  by 

6^  + 1 


where  v is  the  kinematic  viscosity  measured  in  lattice  units. 


(18) 
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3.  Two-dimensional  Homogeneous  Isotropic  Turbulence 

3.1.  INITIAL  AND  BOUNDARY  CONDITIONS 


The  initial  condition  of  the  vorticity  is  randomly  determined  by  satisfying 
the  relation, 


E{k)  = \ | u(ki,k2)  |2  /fc'2  = |/cexp 

" If./  1.1^1  /o  **  \ v / 


(19) 


|fc'-fc|<l/2 


where  u denotes  the  vorticity  in  the  Fourier  space,  k'2  = k'f  + k2,  and 
ki  and  &2  are  the  wave  numbers.  The  periodic  boundary  conditions  are 
imposed  in  the  x and  y directions.  The  computational  domain  is  square, 
(0,0)  < ( x,y ) < (2tt,  2tt). 


3.2.  LOW  REYNOLDS  NUMBER  SIMULATION 

In  order  to  compare  the  results  of  the  LBGK  method  with  those  of  the 
MOL,  numerical  simulation  using  the  square  lattice  model  is  carried  out 
for  a low  Reynolds  number.  The  kinematic  viscosity  is  chosen  as  v = 0.01. 
The  number  of  lattice  nodes  is  129  x 129  (129  lattice  nodes  and  128  lattice 
units  in  one  side).  In  this  case,  the  initial  integral  scale  Reynolds  number 
Rl  corresponds  to  Rl  = 31.9,  which  is  expressed  as  Rl  = 51/W/3.  51  and 
r)  denote  the  total  energy  and  the  enstrophy  dissipation  rate,  which  are 
defined  as 

roc  roc 

51  = / E(k)dk,  r?  = 2i/  / kAE{k)dk  (20) 

Jo  Jo 

Perhaps  the  most  striking  verification  of  the  accuracy  of  the  LBGK  method 
is  found  in  the  direct  comparison  of  contour  plots  of  vorticity  between  the 
LBGK  method  and  higher-order  MOL  at  the  same  physical  time.  Fig.4 
displays  comparison  of  vorticity  contours  at  t — 2.0  between  the  square 
lattice  BGK  method  (solid  line)  and  lOth-order  MOL  (dashed  line).  The 
vorticity  distribution  is  extremely  similar  down  to  detailed  structure  in  the 
two  simulations.  Simulation  by  using  the  triangular  (FHP)  lattice  BGK 
method  is  also  performed  for  the  same  initial  condition  on  lattice  nodes.  In 
Fig.5  comparison  of  vorticity  contour  plots  at  is  shown  between  the  square 
lattice  (solid  line)  and  the  FHP  lattice  (dashed  line).  Once  again,  the  plots 
from  the  two  methods  show  excellent  agreement.  In  order  to  investigate 
the  behavior  of  statistical  quantities,  time  history  of  (a)  the  total  energy 
51,  (b)  the  enstrophy  Q and  (c)  the  enstrophy  dissipation  rate  rj  is  shown  in 
Fig.6,  where  the  enstrophy  Q is  defined  as  Q = J0°°  k2E(k)dk.  Here  51  and 
Q are  inviscid  invariants,  and  therefore  are  monotonically  decreasing  in  this 
dissipative  simulation.  The  solid  line,  dotted  line,  and  dashed  line  indicate 
to  these  quantities  for  the  square  lattice  model,  FHP  lattice  model,  and 
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lOth-order  MOL,  respectively.  What  is  evident  in  Fig.6  is  that  the  LBGK 
method  tracks  the  higher-order  MOL  closely  with  respect  to  evolutions  of 
0,  Q and  r/.  Wave  number  spectra  of  the  total  energy  are  also  compared  in 
Fig.7  which  clearly  displays  the  excellent  match  between  the  three  methods. 
In  the  simulation  of  low  Reynolds  number  case,  no  significant  difference  is 
observed  between  the  square  and  triangular  lattice  models  for  this  resolu- 
tion. 


3.3.  HIGH  REYNOLDS  NUMBER  SIMULATION 

As  a large-scale  direct  numerical  simulation  of  high  Reynolds  number  ho- 
mogeneous isotropic  turbulence,  simulation  for  the  case  with  v = 0.0001 
is  carried  out.  This  corresponds  to  the  initial  integral  scale  Reynolds  num- 
ber Ri  = 25500.  The  number  of  lattice  nodes  is  1025  x 1025.  Fig.8  shows 
comparison  of  vorticity  contour  plots  at  t = 3.0  between  the  square  lat- 
tice BGK  method  and  the  lOth-order  MOL.  Although  slight  difference  in 
vorticity  contours  is  noticeable  at  late  time,  strikingly  similar  features  can 
be  found  for  the  LBGK  simulation,  as  compared  with  the  solutions  by  the 
lOth-order  MOL. 

Time  history  of  (a)  the  total  energy  fl  and  (b)  the  enstrophy  dissipation 
rate  rj  is  shown  in  Fig.9.  Since  Reynolds  number  is  much  higher  than  the 
previous  case,  decrease  in  Q is  much  less  than  observed  in  Fig.6(a).  In 
contrast  to  O and  Q which  are  monotonically  decreasing  in  these  cases,  r\ 
can  be  amplified,  as  much  as  dissipated,  and  not  monotonic  as  shown  in 
Fig.9(b).  Wave  number  spectra  of  k3E(k)  at  t = 3.0  are  compared  in  Fig.10. 
From  this  figure  it  is  seen  that  two  methods  yield  quite  a similar  answer 
in  terms  of  the  statistical  behavior  of  the  flow.  With  the  present  lattice  of 
1025  x 1025  nodes,  the  inertial  range  of  two-dimensional  turbulence  can  be 
resolved.  The  spectrum  shown  in  Fig.10  indicates  that  there  is  a range  of 
wave  number  k < 50  for  which  k3E(k)  is  roughly  constant  so  that  E(k)  is 
proportional  to  k~3. 

Computational  cost  of  the  LBGK  method  for  this  high  Reynolds  number 
simulation  on  SGI  POWER  ONYX  10000  is  compared  in  table  1 with  that 
of  the  lOth-order  MOL.  As  far  as  efficiency  is  concerned,  the  LBGK  method 
requires  less  than  half  CPU  time  per  characteristic  time  of  that  of  the  MOL. 


4.  Three-dimensional  Homogeneous  Isotropic  Turbulence 
4.1.  INITIAL  AND  BOUNDARY  CONDITIONS 


Three-dimensional  decaying  turbulence  is  simulated  with  a random  initial 
condition  having  the  energy  spectrum; 


E(k)  = 16^2 7 


nvlkJaxk^eyiy  -2 {k/kmax)2 


(21) 
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where  we  set  uq  = 1.0  and  the  peak  wave  number  kmax  = 4.75683  for  low 
Reynolds  number  case  and  kmax  = 2.37841  for  medium  Reynolds  num- 
ber one.  The  periodic  boundary  conditions  are  imposed  in  the  x,  y and  z 
directions.  The  computational  domain  is  a cube,  0 < ( x , y,  z)  < 2n. 

4.2.  LOW  AND  MEDIUM  REYNOLDS  NUMBER  SIMULATIONS 

Numerical  simulations  are  carried  out  for  two  cases.  In  the  first  case  the 
kinematic  viscosity  is  chosen  as  v = 0.01189.  The  initial  integral  scale  and 
micro  scale  Reynolds  number  is  30  ~ 45  which  corresponds  to  low  Reynolds 
number  simulation.  In  the  medium  Reynolds  number  case  v — 0.0025  is 
chosen  which  corresponds  the  initial  integral  scale  and  micro  scale  Reynolds 
number=340.  The  number  of  lattice  nodes  is  65  x 65  x 65. 

Contour  surface  of  the  enstrophy  at  t — 0.5  for  the  cubic  LBGK  method 
and  the  MOL  for  v — 0.01189  are  shown  in  Fig.  11  (a)  and  (b)  respectively. 
The  enstrophy  Q is  defined  as 

Q(x,y,z)  = +u>1  +u>2z)  (22) 

Plot  value  is  i\i„t  = 150.  Contour  surface  of  the  enstrophy  computed  by 
two  methods  are  almost  indistinguishable.  Time  history  of  the  total  energy 
and  the  enstrophy  for  medium  Reynolds  number  are  compared  in  Fig.  12(a) 
and  (b).  The  former  is  monotorically  decreasing  while  the  latter  is  ampli- 
fied. Wave  number  dependence  of  energy  spectrum  E(k)  at  t = 0.7  is  also 
compared  in  Fig.13.  From  these  figures  it  is  seen  that  two  methods  yield 
quite  a similar  results  in  terms  of  the  statistical  behavior  of  the  flow.  Com- 
putational cost  of  the  LBGK  method  for  low  Reynolds  number  simulation 
on  a SGI  POWER  ONYX  10000  is  also  compared  in  Table  1 with  that  of 
the  lOth-order  MOL.  As  far  as  efficiency  is  concerned,  the  LBGK  method 
requires  26%  less  CPU  time  than  that  of  the  MOL.  Comparison  between 
the  cubic  LBGK  method  and  the  MOL  shows  that  the  cubic  LBGK  method 
can  be  an  alternative  for  solving  the  Navier-Stokes  equations. 

5.  Conclusions 

Two-  and  three-dimensional  simulations  of  decaying  homogeneous  isotropic 
turbulence  using  the  LBGK  method  have  shown  that  the  method  is  as  ac- 
curate as  the  conventional  method  using  the  same  lattice  size.  The  LBGK 
method  is  able  to  reproduce  the  dynamic  of  decaying  turbulence  and  could 
be  an  alternative  for  solving  the  Navier-Stokes  equations.  Further  investi- 
gation is  needed  on  the  accuracy  and  efficiency  of  cubic  LBGK  model. 
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Figure  1.  Square  lattice  model 


Figure  2.  Triangular  lattice  model 
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Figure  3.  Cubic  lattice  model  countours  between  LBGK(square)  and 

MOL(lOth-order)  at  t = 2.0,  v = 0.01 


LBGK(square) 

LBGK(iriangular) 


Figure  5.  Comparison  of  vorticity 
countours  between  LBGK(square)  and 
LBGK(triangular)  at  t = 2.0,  v — 0.01 


Figure  6(b).  Time  history  of  enstrophy 
Q computed  by  LBGK(square), 
LBGK(triangular)  and  MOL(lOth-order) 
for  v — 0.01 


Figure  6(a).  Time  history  of  total 
energy  fl  computed  by  LBGK(square), 
LBGK(triangular)  and  MOL(lOth-order) 
for  v = 0.01 


Time  !U 


Figure  6(c).  Time  history  of  enstrophy 
disspation  rate  q computed  by 
LBGK(square),  LBGK(triangular)  and 
MOL(lOth-order)  for  v — 0.01 
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Figure  7.  Wave  number  dependence  of 
energy  spectrum  E(k)  computed  by 
LBGK(square),  LBGK(triangular)  and 
MOL(lOth-order)  at  t — 1.0,  v = 0.01 


• i 4 < i h n 

Time  HI 

Figure  9(a).  Time  history  of  total 
energy  Q computed  by  LBGK(square) 
and  MOL(lOth-order)  for  v = 0.001 


Figure  1 0.  Wave  number  dependence  of 
energy  spectrum  kiE(k)  computed  by 
LBGK(square)  and  MOL(lOth-order)  at 
t = 3.0,  v = 0.0001 
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Figure  8.  Comparison  of  vorticity 
countours  between  LBGK(square)  and 
MOL(lOth-order)  at  t — 3.0,  v = 0.0001 


Figure  9(b).  Time  history  of  energy 
dissipation  rate  r]  computed  by 
LBGK(square)  and  MOL(lOt.h-order)  for 
u = 0.001 


Figure  11(a).  Enstorophy  contours  at 
t = 0.5  for  the  cubic  LBGK  method 
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